function MertonModelImpliedVols
clc;clf;close all;

% For the COS method

tau    = [2]';
r      = 0;

% Merton model parameters

S0     = 100;
CP     = 'c';
sigma  = 0.25;
muJ    = 0.0;
sigmaJ = 0.2;
xiP    = 0.1;

% Range of strike prices

K      = linspace(40,180,10)'; 

% Effect of sigmaJ

sigmaJV =[0.0, 0.1, 0.2, 0.3, 0.4];
ivSigmaM = zeros(length(K),length(sigmaJV));
for i=1:length(sigmaJV)
    sigmaJtemp = sigmaJV(i);
    callPrice = MertonCallPrice(CP,S0,K,r,tau,muJ,sigmaJtemp,sigma,xiP);
    for j = 1:length(K)
        strike = K(j);
        call   = callPrice(j);
        ivSigmaM(j,i) = ImpliedVolatility(CP,call,strike,tau,S0,r,0.3)*100;
    end
end
MakeFigure_sigma_J(K, ivSigmaM)

% Effect of xiP

xiPV =[0.0, 0.25, 0.5, 0.75, 1.0];
ivxiPM = zeros(length(K),length(xiPV));
for i=1:length(xiPV)
    xiPtemp = xiPV(i);
    callPrice = MertonCallPrice(CP,S0,K,r,tau,muJ,sigmaJ,sigma,xiPtemp);
    for j = 1:length(K)
        strike = K(j);
        call   = callPrice(j);
        ivxiPM(j,i) = ImpliedVolatility(CP,call,strike,tau,S0,r,0.3)*100;
    end
end
MakeFigure_xiP(K, ivxiPM)

% Effect of muJ

muJV =[-0.25, -0.05, 0.1, 0.25, 0.4];
ivmuJPM = zeros(length(K),length(xiPV));
for i=1:length(xiPV)
    muJtemp = muJV(i);
    callPrice = MertonCallPrice(CP,S0,K,r,tau,muJtemp,sigmaJ,sigma,xiP);
    for j = 1:length(K)
        strike = K(j);
        call   = callPrice(j);
        ivmuJPM(j,i) = ImpliedVolatility(CP,call,strike,tau,S0,r,0.3)*100;
    end
end
MakeFigure_muJ(K, ivmuJPM)

function valueExact = MertonCallPrice(CP,S0,K,r,tau,muJ,sigmaJ,sigma,xiP)
X0  = log(S0);

% Term for E(exp(J)-1)

helpExp = exp(muJ + 0.5 * sigmaJ * sigmaJ) - 1.0;
           
% Analytic representation for Merton's option price

muX     = @(n) X0 + (r - xiP * helpExp - 0.5 * sigma * sigma) * tau + n * muJ;
sigmaX  = @(n) sqrt(sigma * sigma + n * sigmaJ * sigmaJ / tau); 
d1      = @(n) (log(S0./K) + (r - xiP * helpExp - 0.5*sigma * sigma ...
         + sigmaX(n)^2.0) * tau + n * muJ) / (sigmaX(n) * sqrt(tau));
d2      = @(n) d1(n) - sigmaX(n) * sqrt(tau);
value_n = @(n) exp(muX(n) + 0.5*sigmaX(n)^2.0*tau)...
          * normcdf(d1(n)) - K .* normcdf(d2(n));
      
% Option value calculation, it is an infinite sum but we truncate at 20 terms

valueExact = value_n(0.0);
kidx = 1:20;
for k = kidx
    valueExact = valueExact+  (xiP * tau).^k * value_n(k) / factorial(k);
end
valueExact = valueExact * exp(-r*tau) * exp(-xiP * tau);

if CP == 'p'
    valueExact = valueExact - S0 + K * exp(-r*tau);
end

% Closed-form expression of European call/put option with Black-Scholes formula

function value=BS_Call_Option_Price(CP,S_0,K,sigma,tau,r)

% Black-Scholes call option price

d1    = (log(S_0 ./ K) + (r + 0.5 * sigma^2) * tau) / (sigma * sqrt(tau));
d2    = d1 - sigma * sqrt(tau);
if lower(CP) == 'c' || lower(CP) == 1
    value =normcdf(d1) * S_0 - normcdf(d2) .* K * exp(-r * tau);
elseif lower(CP) == 'p' || lower(CP) == -1
    value =normcdf(-d2) .* K*exp(-r*tau) - normcdf(-d1)*S_0;
end

function impliedVol = ImpliedVolatility(CP,marketPrice,K,T,S_0,r,initialVol)
    func = @(sigma) (BS_Call_Option_Price(CP,S_0,K,sigma,T,r) - marketPrice).^1.0;
    impliedVol = fzero(func,initialVol);
    
function MakeFigure_sigma_J(X1, YMatrix1)

%CREATEFIGURE(X1,YMATRIX1)
%  X1:  vector of x data
%  YMATRIX1:  matrix of y data

%  Auto-generated by MATLAB on 16-Jan-2012 15:26:40

% Create figure

figure1 = figure('InvertHardcopy','off',...
    'Colormap',[0.061875 0.061875 0.061875;0.06875 0.06875 0.06875;0.075625 0.075625 0.075625;0.0825 0.0825 0.0825;0.089375 0.089375 0.089375;0.09625 0.09625 0.09625;0.103125 0.103125 0.103125;0.11 0.11 0.11;0.146875 0.146875 0.146875;0.18375 0.18375 0.18375;0.220625 0.220625 0.220625;0.2575 0.2575 0.2575;0.294375 0.294375 0.294375;0.33125 0.33125 0.33125;0.368125 0.368125 0.368125;0.405 0.405 0.405;0.441875 0.441875 0.441875;0.47875 0.47875 0.47875;0.515625 0.515625 0.515625;0.5525 0.5525 0.5525;0.589375 0.589375 0.589375;0.62625 0.62625 0.62625;0.663125 0.663125 0.663125;0.7 0.7 0.7;0.711875 0.711875 0.711875;0.72375 0.72375 0.72375;0.735625 0.735625 0.735625;0.7475 0.7475 0.7475;0.759375 0.759375 0.759375;0.77125 0.77125 0.77125;0.783125 0.783125 0.783125;0.795 0.795 0.795;0.806875 0.806875 0.806875;0.81875 0.81875 0.81875;0.830625 0.830625 0.830625;0.8425 0.8425 0.8425;0.854375 0.854375 0.854375;0.86625 0.86625 0.86625;0.878125 0.878125 0.878125;0.89 0.89 0.89;0.853125 0.853125 0.853125;0.81625 0.81625 0.81625;0.779375 0.779375 0.779375;0.7425 0.7425 0.7425;0.705625 0.705625 0.705625;0.66875 0.66875 0.66875;0.631875 0.631875 0.631875;0.595 0.595 0.595;0.558125 0.558125 0.558125;0.52125 0.52125 0.52125;0.484375 0.484375 0.484375;0.4475 0.4475 0.4475;0.410625 0.410625 0.410625;0.37375 0.37375 0.37375;0.336875 0.336875 0.336875;0.3 0.3 0.3;0.28125 0.28125 0.28125;0.2625 0.2625 0.2625;0.24375 0.24375 0.24375;0.225 0.225 0.225;0.20625 0.20625 0.20625;0.1875 0.1875 0.1875;0.16875 0.16875 0.16875;0.15 0.15 0.15],...
    'Color',[1 1 1]);

% Create axes

%axes1 = axes('Parent',figure1,'Color',[1 1 1]);
axes1 = axes('Parent',figure1);
grid on

% Uncomment the following line to preserve the X-limits of the axes
% xlim(axes1,[45 160]);
% Uncomment the following line to preserve the Y-limits of the axes
% ylim(axes1,[19 26]);
% Uncomment the following line to preserve the Z-limits of the axes
% zlim(axes1,[-1 1]);

box(axes1,'on');
hold(axes1,'all');

% Create multiple lines using matrix input to plot
% plot1 = plot(X1,YMatrix1,'Parent',axes1,'MarkerEdgeColor',[0 0 0],...
%     'LineWidth',1,...
%     'Color',[0 0 0]);

plot1 = plot(X1,YMatrix1,'Parent',axes1,...
    'LineWidth',1);
set(plot1(1),'Marker','diamond','DisplayName','\sigma_J=0.0');
set(plot1(2),'Marker','square','LineStyle','-.',...
    'DisplayName','\sigma_J=0.1');
set(plot1(3),'Marker','o','LineStyle','-.','DisplayName','\sigma_J=0.2');
set(plot1(4),'DisplayName','\sigma_J=0.3');
set(plot1(5),'Marker','v','DisplayName','\sigma_J=0.4');

% Create xlabel

xlabel({'K'});

% Create ylabel

ylabel({'implied volatility [%]'});

% Create title

title({'Effect of \sigma_J on implied volatility'});

% Create legend

legend1 = legend(axes1,'show');
set(legend1,'Color',[1 1 1]);

function MakeFigure_xiP(X1, YMatrix1)

%CREATEFIGURE(X1,YMATRIX1)
%  X1:  vector of x data
%  YMATRIX1:  matrix of y data

%  Auto-generated by MATLAB on 16-Jan-2012 15:26:40

% Create figure

figure1 = figure('InvertHardcopy','off',...
    'Colormap',[0.061875 0.061875 0.061875;0.06875 0.06875 0.06875;0.075625 0.075625 0.075625;0.0825 0.0825 0.0825;0.089375 0.089375 0.089375;0.09625 0.09625 0.09625;0.103125 0.103125 0.103125;0.11 0.11 0.11;0.146875 0.146875 0.146875;0.18375 0.18375 0.18375;0.220625 0.220625 0.220625;0.2575 0.2575 0.2575;0.294375 0.294375 0.294375;0.33125 0.33125 0.33125;0.368125 0.368125 0.368125;0.405 0.405 0.405;0.441875 0.441875 0.441875;0.47875 0.47875 0.47875;0.515625 0.515625 0.515625;0.5525 0.5525 0.5525;0.589375 0.589375 0.589375;0.62625 0.62625 0.62625;0.663125 0.663125 0.663125;0.7 0.7 0.7;0.711875 0.711875 0.711875;0.72375 0.72375 0.72375;0.735625 0.735625 0.735625;0.7475 0.7475 0.7475;0.759375 0.759375 0.759375;0.77125 0.77125 0.77125;0.783125 0.783125 0.783125;0.795 0.795 0.795;0.806875 0.806875 0.806875;0.81875 0.81875 0.81875;0.830625 0.830625 0.830625;0.8425 0.8425 0.8425;0.854375 0.854375 0.854375;0.86625 0.86625 0.86625;0.878125 0.878125 0.878125;0.89 0.89 0.89;0.853125 0.853125 0.853125;0.81625 0.81625 0.81625;0.779375 0.779375 0.779375;0.7425 0.7425 0.7425;0.705625 0.705625 0.705625;0.66875 0.66875 0.66875;0.631875 0.631875 0.631875;0.595 0.595 0.595;0.558125 0.558125 0.558125;0.52125 0.52125 0.52125;0.484375 0.484375 0.484375;0.4475 0.4475 0.4475;0.410625 0.410625 0.410625;0.37375 0.37375 0.37375;0.336875 0.336875 0.336875;0.3 0.3 0.3;0.28125 0.28125 0.28125;0.2625 0.2625 0.2625;0.24375 0.24375 0.24375;0.225 0.225 0.225;0.20625 0.20625 0.20625;0.1875 0.1875 0.1875;0.16875 0.16875 0.16875;0.15 0.15 0.15],...
    'Color',[1 1 1]);

% Create axes

%axes1 = axes('Parent',figure1,'Color',[1 1 1]);
axes1 = axes('Parent',figure1);
grid on

% Uncomment the following line to preserve the X-limits of the axes
% xlim(axes1,[45 160]);
% Uncomment the following line to preserve the Y-limits of the axes
% ylim(axes1,[19 26]);
% Uncomment the following line to preserve the Z-limits of the axes
% zlim(axes1,[-1 1]);

box(axes1,'on');
hold(axes1,'all');

% Create multiple lines using matrix input to plot
% plot1 = plot(X1,YMatrix1,'Parent',axes1,'MarkerEdgeColor',[0 0 0],...
%     'LineWidth',1,...
%     'Color',[0 0 0]);

plot1 = plot(X1,YMatrix1,'Parent',axes1,...
    'LineWidth',1);
set(plot1(1),'Marker','diamond','DisplayName','\xi_P=0.0');
set(plot1(2),'Marker','square','LineStyle','-.',...
    'DisplayName','\xi_P=0.25');
set(plot1(3),'Marker','o','LineStyle','-.','DisplayName','\xi_P=0.5');
set(plot1(4),'DisplayName','\xi_P=0.75');
set(plot1(5),'Marker','v','DisplayName','\xi_P=1.0');

% Create xlabel

xlabel({'K'});

% Create ylabel

ylabel({'implied volatility [%]'});

% Create title

title({'Effect of \xi_P on implied volatility'});

% Create legend

legend1 = legend(axes1,'show');
set(legend1,'Color',[1 1 1]);

function MakeFigure_muJ(X1, YMatrix1)

%CREATEFIGURE(X1,YMATRIX1)
%  X1:  vector of x data
%  YMATRIX1:  matrix of y data

%  Auto-generated by MATLAB on 16-Jan-2012 15:26:40

% Create figure

figure1 = figure('InvertHardcopy','off',...
    'Colormap',[0.061875 0.061875 0.061875;0.06875 0.06875 0.06875;0.075625 0.075625 0.075625;0.0825 0.0825 0.0825;0.089375 0.089375 0.089375;0.09625 0.09625 0.09625;0.103125 0.103125 0.103125;0.11 0.11 0.11;0.146875 0.146875 0.146875;0.18375 0.18375 0.18375;0.220625 0.220625 0.220625;0.2575 0.2575 0.2575;0.294375 0.294375 0.294375;0.33125 0.33125 0.33125;0.368125 0.368125 0.368125;0.405 0.405 0.405;0.441875 0.441875 0.441875;0.47875 0.47875 0.47875;0.515625 0.515625 0.515625;0.5525 0.5525 0.5525;0.589375 0.589375 0.589375;0.62625 0.62625 0.62625;0.663125 0.663125 0.663125;0.7 0.7 0.7;0.711875 0.711875 0.711875;0.72375 0.72375 0.72375;0.735625 0.735625 0.735625;0.7475 0.7475 0.7475;0.759375 0.759375 0.759375;0.77125 0.77125 0.77125;0.783125 0.783125 0.783125;0.795 0.795 0.795;0.806875 0.806875 0.806875;0.81875 0.81875 0.81875;0.830625 0.830625 0.830625;0.8425 0.8425 0.8425;0.854375 0.854375 0.854375;0.86625 0.86625 0.86625;0.878125 0.878125 0.878125;0.89 0.89 0.89;0.853125 0.853125 0.853125;0.81625 0.81625 0.81625;0.779375 0.779375 0.779375;0.7425 0.7425 0.7425;0.705625 0.705625 0.705625;0.66875 0.66875 0.66875;0.631875 0.631875 0.631875;0.595 0.595 0.595;0.558125 0.558125 0.558125;0.52125 0.52125 0.52125;0.484375 0.484375 0.484375;0.4475 0.4475 0.4475;0.410625 0.410625 0.410625;0.37375 0.37375 0.37375;0.336875 0.336875 0.336875;0.3 0.3 0.3;0.28125 0.28125 0.28125;0.2625 0.2625 0.2625;0.24375 0.24375 0.24375;0.225 0.225 0.225;0.20625 0.20625 0.20625;0.1875 0.1875 0.1875;0.16875 0.16875 0.16875;0.15 0.15 0.15],...
    'Color',[1 1 1]);

% Create axes

%axes1 = axes('Parent',figure1,'Color',[1 1 1]);
axes1 = axes('Parent',figure1);
grid on

% Uncomment the following line to preserve the X-limits of the axes
% xlim(axes1,[45 160]);
% Uncomment the following line to preserve the Y-limits of the axes
% ylim(axes1,[19 26]);
% Uncomment the following line to preserve the Z-limits of the axes
% zlim(axes1,[-1 1]);

box(axes1,'on');
hold(axes1,'all');

% Create multiple lines using matrix input to plot
% plot1 = plot(X1,YMatrix1,'Parent',axes1,'MarkerEdgeColor',[0 0 0],...
%     'LineWidth',1,...
%     'Color',[0 0 0]);

plot1 = plot(X1,YMatrix1,'Parent',axes1,...
    'LineWidth',1);
set(plot1(1),'Marker','diamond','DisplayName','\mu_J=-0.25');
set(plot1(2),'Marker','square','LineStyle','-.',...
    'DisplayName','\mu_J=-0.05');
set(plot1(3),'Marker','o','LineStyle','-.','DisplayName','\mu_J=0.1');
set(plot1(4),'DisplayName','\mu_J=0.25');
set(plot1(5),'Marker','v','DisplayName','\mu_J=0.4');

% Create xlabel

xlabel({'K'});

% Create ylabel

ylabel({'implied volatility [%]'});

% Create title

title({'Effect of \mu_J on implied volatility'});

% Create legend

legend1 = legend(axes1,'show');
set(legend1,'Color',[1 1 1]);
